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1.  Abstract 

Our  goal  was  to  develop  novel  multiparticle  computational  tools  for  materials 
science  based  on  representing  operators  and  functions  of  many  variables  as  short 
sums  of  separable  functions. 

We  have  extended  multiresolution  separated  representation  of  free-space  Green’s 
functions  to  those  periodic  or  satisfying  boundary  conditions,  making  fast  algo¬ 
rithms  available  for  solving  integral  equations. 

The  problem  of  incorporating  inter-electron  cusps  within  two-particle  methods 
led  us  to  the  problem  of  approximating  multivariable  functions  by  sums  of  prod¬ 
ucts  of  Gaussians.  We  have  developed  and  tested  a  new  (suboptimal)  algorithm, 
sufficient  for  many  practical  purposes.  Using  this  algorithm  we  constructed  novel 
separated  representations  of  non-convolutional  Green’s  functions  and  spectral  pro¬ 
jectors  for  operators  with  potentials  from  the  Rollnik  class.  Such  potentials  include 
all  physically  significant  potentials  considered  within  a  finite  domain.  For  optimal 
representations,  we  have  developed  an  algorithm  in  the  case  of  two  variables;  ex¬ 
tensions  to  higher  dimensions  need  to  be  worked  on  further.  We  started  work  on 
moving  preliminary  implementations  of  these  algoritrhms  into  the  Python  environ¬ 
ment  where  they  can  be  used  for  practical  computations. 

We  have  developed  the  algorithmic  framework  for  a  size-extensive/consistent 
method  for  the  multiparticle  Schrodinger  equation.  The  central  computation  in¬ 
volved  has  revealed  a  center-of-mass  principle  for  electrons,  which  gives  hope  for  the 
construction  of  an  analogue  of  the  Fast  Multipole  Method  for  quantum  mechanical 
systems. 


2.  Tasks 

Our  work  addresses  three  tasks  in  using  separated  representations  [11,  9]  for  com¬ 
putational  materials  science.  In  our  proposal  we  organized  these  tasks  by  increasing 
degree  of  difficulty.  They  are: 

1.  Incorporating  multiresolution  separated  representations  into  existing  methods 
in  materials  science  (where  the  underlying  dimension  is  three)  and,  thus,  sim¬ 
plifying  computations  involving  Green’s  functions  with  boundary  conditions 
(or  periodic) 

2.  “Upgrading”  computational  chemistry  from  one-particle  theories  to  using  two- 
particle  theories  (with  underlying  dimension  six). 
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3.  Developing  methods  for  solving  the  multiparticle  Schrddinger  equation  of 
quantum  mechanics,  where  the  dimension  grows  linearly  with  the  number 
of  particles.  In  particular,  our  goal  is  finding  computational  methods  that 
incorporate  algorithmic  size-extensivity. 

3.  Results  for  Task  1:  Multiresolution  separated  representations  of 
Green’s  functions  with  boundary  conditions 

Our  goal  has  been  to  obtain  multiresolution  representations  of  lattice  sums  for 
Green’s  functions  (as  well  as  potentials)  in  a  form  that  facilitates  solving  inte¬ 
gral  equations.  We  note  that  computations  involving  multidimensional  free  space 
Green’s  functions  are  greatly  simplified  by  using  separated  representations  [11,  9, 
28,  30,  50,  51,  22,  8].  Using  lattice  sums,  we  have  constructed  separated  representa¬ 
tions  of  Green’s  functions  with  periodic,  Dirichlet  or  Neumann  boundary  conditions. 
Such  lattice  sums  are  not  absolutely  convergent  and  we  have  developed  an  explicit 
multiresolution  interpretation  of  these  sums. 

Our  approach  is  based  on  representing  spherically  symmetric  functions  by  sums 
of  products  of  Gaussians.  Such  approximations  are  obtained  by  discretizing  inte¬ 
grals  which  are  similar,  or  even  identical,  to  those  used  in  the  Ewald  summation. 
However,  even  though  conceptually  our  construction  has  strong  similarities  with 
the  Ewald  summation,  it  differs  in  important  details.  We  compute  separated  rep¬ 
resentations  of  lattice  sums  using  only  one  dimensional  integrals  or  sums.  The 
resulting  approximation  is  a  combinations  of  three  terms,  a  smooth  periodic  term 
with  a  short  Fourier  expansion,  and  two  terms,  a  singular  and  non-singular,  the 
sum  of  which  is  near  zero  outside  the  ball  inscribed  into  a  unit  cell  (for  any  finite 
but  arbitrary  accuracy) .  The  number  of  terms  needed  in  the  final  representation  is 
(nearly)  optimal. 

3.1.  Technical  Introduction:  separated  representations  of  lattice  sums. 

A  practical  way  of  computing  Green’s  functions  and  periodic  potentials  via  lattice 
sums  has  been  of  interest  for  many  years.  An  early  seminal  contribution  was  made 
by  P.  Ewald  [21],  although  the  history  of  lattice  sums  starts  earlier  and  we  refer  to 
[25]  for  a  historical  overview  and  results  prior  to  1980.  A  more  recent  fundamental 
advance  was  made  in  [32],  where  a  surprisingly  simple  integral  representations  of 
harmonic  lattice  sums  has  been  derived. 

Our  goal  is  to  compute  representations  of  lattice  sums  for  Green’s  functions 
and  potentials  in  a  form  that  facilitates  solving  integral  equations.  We  note  that 
computations  involving  multidimensional  free  space  Green’s  functions  are  greatly 
simplified  by  using  separated  representations  [9,  10,  28,  30,  50,  51,  22,  8].  Using 
the  lattice  sums,  we  construct  separated  representations  of  Green’s  functions  with 
periodic  or  zero  boundary  conditions.  Such  lattice  sums  are  not  absolutely  conver¬ 
gent  and  we  develop  an  explicit  multiresolution  interpretation  of  these  sums.  We 
note  that  we  selected  the  Green’s  functions  for  the  Poisson  and  the  bound  state 
Helmholtz  operators  as  examples  since  the  same  method  applies  equally  well  to 
other  non-oscillatory  Green’s  functions. 

For  spherically  symmetric  singular  potentials  it  has  been  traditional  to  rely  on 
spherical  harmonics  as  a  tool  for  computing  lattice  sums.  The  essential  difficulty 
in  such  approach  is  to  combine  together  local  spherical  symmetry  and  summation 
over  directional  shifts.  In  our  approach  we  avoid  this  difficulty  by  representing 
spherically  symmetric  functions  in  the  Cartesian  system  of  coordinates  by  using 
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sums  of  products  of  Gaussians.  Such  separated  representations  are  constructed  to 
be  accurate  over  a  remarkably  wide  range  of  spatial  scales  since  the  number  of  terms 
grows  only  as  0(—  log(<5)),  where  6  is  the  distance  away  from  the  singularity.  The 
summation  over  directional  shifts  then  proceeds  in  the  term  by  term  fashion  and 
involves  sums  only  over  one-dimensional  lattices. 

3.2.  Preliminary  Considerations. 


3.2.1.  Separated  representations  for  Poisson-type  kernels.  Let  us  construct  a  sepa¬ 
rated  approximation  of  the  function  l/r“,  where  r  =  ||x||,  x  e  IR.3  using  a  collection 
of  Gaussians.  The  approximation  is  obtained  by  first  discretizing  the  integral 


(3.1) 


0-re  -\-ots 


r(a/2)  J_c 


ds . 


For  a  =  1  it  is  the  same  integral  as  used  in  the  Ewald  summation  (up  to  a  change 
of  variables,  see  e.g.,  [25]).  We  have 


Proposition  3.1.  For  any  a  >  0,  0  <  S  <  1,  and  0  <  e  <  min  {|,  — there  exist 
positive  numbers  pm  and  wm  such  that 


(3.2) 

where 


M 

E 

771=1 


wme 


—Pr, 


e 

-  j.OL  1 


(3.3)  M  =  log e  1  [c0  +  ci  log  e  1+C2log<5  1], 

where  c\,  C2  and  C3  are  constants  that  only  depend  on  a.  For  fixed  power  a  and 
accuracy  e,  we  have  M  =  O(log5_1). 


A  proof  of  Proposition  3.1  can  be  found  in  [13]. 

Using  r  =  ||x||,  where  x  =  (x\,  X2,  £3),  and  a  =  1  in  (3.2),  we  arrive  at  a 
separated  representation  for  the  Poisson  kernel.  Although  in  this  paper  we  compute 
the  lattice  sums  corresponding  to  the  Poisson  kernel,  the  same  approach  will  work 
for  any  a  >  0  as  well  as  other  spherically  symmetric  potentials,  e.g.  Yukawa 
potential  e_/ir/r. 

As  in  [28,  30],  the  approximation  in  (3.2)  is  obtained  using  trapezoidal  rule.  First, 
we  discretize  the  integral  (3.1),  namely,  set  pm  =  e2sm  and  wm  =  2A s  eaSm /T(a/2), 
where  sm  =  so  +  (to  —  l)As,  to  =  1 . . . ,  M.  For  a  given  accuracy  e  and  range 
0  <  5  <  r  <  1,  we  select  so  and  sm  =  sq  +  ( M  —  l)As,  the  end  points  of  the 
interval  of  integration  replacing  the  real  line  in  (3.1),  so  that  at  these  points  the 
function  f(s)  =  e~r  e  +as  and  a  sufficient  number  of  its  derivatives  are  close  to 
zero  to  within  the  desired  accuracy.  We  also  select  M,  the  number  of  points  in  the 
quadrature,  so  that  the  accuracy  requirement  is  satisfied. 


3.2.2.  The  cross- correlation  functions.  We  use  the  scaling  functions  of  the  multi- 
wavelet  bases  developed  in  [2].  For  a  brief  review  of  the  multiwavelet  bases  see 
also  [3].  For  convolution  operators  we  only  need  to  compute  integrals  with  the 
cross-correlation  functions  of  the  scaling  functions,  namely, 


0  <  2;  <  1, 
— 1  <  2;  <  0, 
0,  1  <  |2;|, 


(3.4) 


<E hv  (x) 
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where  i,i'  —  0, . . .  ,  m  —  1,  m  is  the  order  of  the  basis,  and 

pl  —  X  rO 

(3-5)  $u,{x)  =  /  <j>i(x  +  y)<j>i'{y)dy,  $^,(ar)  =  /  4>i{x  +  y)  <j>i>{y)dy . 

Jo  j  -x 

The  scaling  functions  fa  are  the  normalized  Legendre  polynomials  on  the  interval 
[0,1], 

A  (T\  _  I  v/2lzTlPi(2x  -  1),  X  e  [0, 1] 

M  j"l  0,  x$[0,l]  ’ 

where  Pi  are  the  Legendre  polynomials  on  [—1, 1],  This  implies  that  the  functions 
are  piecewise  polynomials  of  degree  i  +  i'  +  1  with  the  support  in  [—1,1]. 

Proposition  3.2. 

1.  Transposition  of  indices:  &a>(x)  =  (— l)I+i  $i>i(x), 

2.  Relations  between  $+  and  . x)  =  (— 1)*+®  for  0  <  x  <  1, 

3.  Values  at  zero:  3>u>  (0)  =  0  if  i  ^  i' ,  and  $^(0)  =  1  for  i  =  0, . . . ,  m  —  I* 

4.  Upper  bound:  maxx€r_lii]  |$jj/(a;)|  <  1  for  i,  i'  =  0, . . . ,  m  —  1, 

5.  Connection  with  the  Gegenbauer  polynomials: 

^0{x)  =  i  c[~1/2\ 2x  -  1)  +  \  and  $+(*)  =  \  y/W+lC^/2)  (2x  -  1),  for 
l  —  1,2, ,  where  is  the  Gegenbauer  polynomial, 

6.  Linear  expansion:  ( a;)  are  linear  combinations 


(3.6) 

^'0*0=  X!  4  *«(*)> 

where 

(3.7) 

4'=/  ®tAx)®io(x)0- x2)~%dx 

Jo 

7.  Vanishing  moments:  we  have  $oo(a;)  dx  =  1  and  xk^u>(x)  dx  =  0  for 
i  +  i'  >  1  and  D  <k  <  i  +  i'  —  1. 

3.3.  Multiresolution  representation  of  lattice  sums.  Since  the  sum  in  (3.14) 
is  not  absolutely  convergent,  we  need  to  provide  its  interpretation.  Our  derivation 
in  previous  sections  used  the  usual  justification  for  Ewald  summation,  see  e.g.  [25]. 
In  this  section  we  interpret  the  divergent  sum  in  (3.14)  using  multiresolution  repre¬ 
sentation  of  the  Poisson  kernel  as  a  starting  point.  This  allows  us  to  obtain  simple 
and  transparent  conditions  for  the  shape  of  the  domain  over  which  the  summation 
is  performed  as  well  as  to  observe  that  the  complexity  of  the  algorithm  to  apply  the 
periodic  Green’s  function  is  the  same  as  that  for  the  Green’s  function  in  R3  where 
the  differences  between  the  two  occurs  only  on  coarse  scales. 

Let  us  begin  with  the  multiresolution  representation  of  the  Poisson  kernel  in 
R3.  We  use  multiwavelet  bases  [2]  as  we  find  them  more  appropriate  for  numerical 
applications  [3] ,  [30]  although  other  wavelet  bases  may  be  used  in  this  construction 
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as  well.  The  multiwavelet  bases  are  constructed  in  L2{B),  B  =  [— 1/2,  l/2]3and 
then  extended  to  form  a  basis  in  T2(R3)  by  replicating  the  construction  for  each 
cube  shifted  by  n  £  Z3  (see  e.g.  [3]).  Let  Vj  €  T2(R3)  denote  subspaces  of  the 
multiresolution  analysis  corresponding  to  a  multiwavelet  basis  of  order  m,  where  m 
is  the  number  of  nodes  in  the  Gaussian  quadrature  (in  each  direction).  We  then 
have 


V0  C  Vi  C  V2  C  ■  ■  ■  C  Vj  ■  ■  ■ 

or,  for  each  Vn, 

Vn  =  V0  ©  W0  ©  Wx  ©  •  •  •  ©  W„_! 

where  Vj  ®  Wj  =  Vj+i-  The  subspace  Vo  is  spanned  by  products  of  orthogonal 
polynomials  up  to  degree  m  —  1  in  each  variable  localized  within  cubes  shifted  by 

n  e  z3. 

Let  Pj  and  Qj  denote  the  orthogonal  projectors  on  spaces  Vj  and  Wj,  respec¬ 
tively.  Given  an  operator  T ,  we  obtain  its  telescopic  series  as 

T  =  T0  +  (Ti  -  T0)  +  (T2  -  Ti)  +  . . . , 

where  Tj  =  PjTPj  is  the  projection  of  the  operator  T  on  the  subspace  Vj.  It  is 
easy  to  see  that 

Tj+i  —  Tj  =  QjTQj  +  QjTPj  +  PjTQj , 

is  the  element  of  the  non-standard  form  of  the  operator  T  on  the  scale  j  [6] .  It  was 
shown  in  [6]  that  the  vanishing  moments  of  the  wavelet  basis  (this  property  remains 
unchanged  for  multiwavelets)  yield  a  rapid  decay  of  the  coefficients  of  Tj+\—Tj  away 
from  the  singularity  of  the  kernel.  The  rate  of  decay  is  controlled  by  the  number 
of  the  vanishing  moments  which  we  choose  depending  on  the  accuracy  requirement 
as  will  be  explained  below.  By  choosing  m  >  2,  we  assure  the  rate  of  decay  of  at 
least  l/rm+1  for  all  terms  in  the  telescopic  series  except  for  T0.  If  we  now  sum  over 
the  periodic  lattice,  separately  on  each  scale,  the  sum  is  rapidly  convergent  for  all 
terms  except  that  for  To-  In  fact,  by  selecting  the  number  of  vanishing  moments  of 
the  basis,  we  control  the  rate  of  decay  and  may  choose  that  number  so  that  for  a 
fixed  but  arbitrary  precision  the  contribution  of  all  outside  terms  is  negligible. 

What  remains  to  be  shown  is  how  to  perform  the  summation  on  Vq,  making 
it  the  key  to  the  multiresolution  definition  of  the  periodic  Green’s  function.  The 
projection  To  of  the  Poisson  kernel  is  given  by  the  integrals 

(3.8)  tu',jj',kk'  =  /  1 1  1 1  l)  $»'0e2)  Qkk'ixz)  dx, 

JB  llx  +  nl! 

where  n  G  Z3,  B  =  [—1/2, 1/2]3,  x  =  (x\,  X2,  £3),  j1,  k,  k!  =  0, . . .  ,  m  —  1  and 

m  is  the  order  of  the  basis  (the  number  of  nodes  in  the  Gaussian  quadrature) .  We 
now  need  to  interpret 

(3.9) 

/  TT— Jl$ii'(xi)$jj'{x2)$kk'{X3)dxL, 
neZ3  neZ3  Jb  11  11 

integrals  with  cross-correlation  functions  rather  than  the  sum  in  (3.14). 
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Let  us  consider  t",  •  kk,  as  a  function  of  n  G  Z3.  Using 


(3.10) 


1 

I  lx  +  n 


!_  _  1  1  |x|  |2  3  (x  ■  n)2  1 

n||  1 1 n| | 3  2  | |n| |3  2  ||n||5  ||n||4 


we  observe  that  the  coefficient  t^0Q  oo  cannot  be  given  any  meaning  through  the 
sum  in  (3.9).  However,  if  we  apply  the  resulting  operator  to  periodic  functions  with 
zero  mean,  then  the  coefficient  tgo'oo  oo  is  n°t  needed  and  we  simply  set  00  =  0 
rather  than  use  (3.9). 

For  i  +  i'  >  1  all  functions  $ij/  (see  Section  3.2.2)  have  vanishing  moments, 
namely, 

(3.11)  J  &u’(x)  xm  dx  =  0,  0<m<i  +  i'  —  1. 


Using  (3.11)  and  (3.10),  integrals  (3.9)  involving  functions  any  function  $ri/  with 
indices  i  +  i'  >  3  result  in  a  rapid  decay  of  the  coefficients  jj,  kk,  =  G(jjT^pr) 
and,  therefore,  the  corresponding  sum  in  (3.9)  will  converge  absolutely. 

Let  us  now  examine  integrals  (3.9)  involving  function  $jj/  with  indices  1  <  i+i'  < 
2.  We  find  that  not  all  slowly  decaying  terms  in  (3.10)  vanish.  For  example,  using 
(3.9)  and  (3.10),  and  the  fact  that  $0i  is  an  odd  function,  we  have 

(3.12)  doi,oo,oo(n)  =  rrri3  (  [  x<&m{x)dx)((  $00(x)dx)2  +  ©(t— 1 -r^), 

llnll  J- i  J-i  llnH 


and,  using  the  fact  that  $n  is  an  even  function, 
(3.13) 

Jii,oo,oo(n)  =  - — -  (  /  a;2  $n (x)dx) 

2  n  5  7-i 


$00  (x)  dx) 2  +  0( 


)• 


In  all  of  these  cases  (conditional)  summation  over  symmetric  ranges  of  indices  can¬ 
cels  the  first  term  and  yields  convergent  sums  for  the  coefficients  kk,,  thus 

completing  the  multiresolution  definition  of  the  periodic  Green’s  function. 


3.4.  Lattice  sums  for  3D  cubic  lattice.  Let  us  now  describe  how  to  use  (3.2) 
in  a  more  traditional  approach  to  lattice  sums.  Let  us  consider  the  cubic  lattice 
(with  some  modifications  the  results  can  be  extended  to  a  general  lattice).  Starting 
with  the  free  space  Poisson  kernel  (or,  alternatively,  the  Coulomb  potential),  let  us 
compute  for  x  €  B, 


(3.14) 


G'o(x)  —  ^2 


■n 


1 


+  £'l 

nez3 


x  +  n 


where  the  unit  cell  B  =  [—1/2, 1  /2] 3 .  Formally  Go  is  the  Green’s  function  in  B, 

— V2Gq(x)  =  47n5(x), 


with  the  periodic  boundary  conditions.  Since  the  sum  (3.14)  is  not  absolutely 
convergent,  summation  in  (3.14)  needs  a  separate  description.  If  we  apply  the 
Green’s  function  Go  only  to  periodic  functions  integrating  to  zero  in  B ,  fB  f(x)dx  = 
0,  then  (3.14)  is  a  conditionally  convergent  sum.  Let  us  first  assume  that  the 
convergence  issues  are  resolved  via  a  procedure  described  in  e.g.  [25]  or  via  a 
multiresolution  justification  of  such  summation  in  Section  3.3. 
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Figure  3.1.  Relative  error  (log10)  of  approximating  the  Poisson 
kernel  in  (3.15),  where  e  =  10— 9 ,  1/2  <  ||x||  <  1015  . 


Let  r o  =  1/2  be  the  radius  of  the  ball  inscribed  into  B  .  We  rescale  (3.2)  to 
approximate  1/r  outside  the  ball  B ,  in  the  region  r  €  [ro,I?],  where  R  is  a  large 
number.  We  obtain  for  1/2  <r<  1/(25), 


(3.15) 


m= 1 


where  tm  =  4 pm52  and  pm  =  2 wm5.  The  error  e(r)  =  r|£  -  X)m=i  Pme~tmr2  |  is 
illustrated  in  Figure  3.1.  In  this  case  M  =  270,  tm  =  e2Tm  and  pm  =  &  •  eTm,  where 
t0  =  -59.004022799786696,  A  =  0.22676579925650819,  cr  =  0.25587780369080371, 
and  Tm  =  to  +  A(m  —  1)  with  m  =  1, . . .  M.  As  it  turns  out,  we  need  only  some 
of  these  terms  and  we  let  the  algorithm  select  the  necessary  ones.  We  note  that  in 
this  example  the  largest  exponent  in  (3.15)  is  «  54. 

Assuming  that  we  assign  the  lattice  sum  S'  a  finite  value  following,  for  example, 
the  recipe  in  [25],  we  have  formally 


(3.16) 


£' 


x  +  n 


m=  1 


£'f 


-tm  x+n 


E'i 


x  +  n 


Using  (3.14)  and  (3.16),  we  obtain  an  approximation  to  G o  as 


M 


Pm& 


G“ppr(x) 


X 


+  CW, 
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where  G^er  is  the  periodic  component, 


M 

(3.17)  C(x)  =  E  E  e-‘-||x+n|12. 

m—  1  ngZ3 


By  construction,  for  x  £  B  the  combination 

M 

(3.18)  n-rr  -  E  Prne~tmlM\ 

II  m—  1 

is  less  than  e  outside  of  the  ball  of  radius  1/2.  Evaluating  the  sum  in  (3.18)  only  for 
x  £  73,  we  significantly  reduce  the  necessary  number  of  terms,  since  within  accuracy 
e  contributions  of  many  terms  in  (3.15)  are  accounted  for  by  a  constant. 

Next  let  us  compute  the  Fourier  coefficients  of  the  periodic  function  Gger.  We 
have  for  p  £  Z3, 


r  M 

9p=  ( y ^  Pm  e 

B  m=l  nGZ3 

and,  therefore, 

(3.19) 


M 

pdx  =  e  p„ 


-tm  x  ^—2nix- 


pdx, 


m=  1 


M 

m=l 


The  coefficient  go  is  not  well-defined  by  (3.19),  since  as  we  improve  the  approxima¬ 
tion  (3.15),  the  expression  (3.19)  diverges.  However,  we  only  apply  Go  to  periodic 
functions  with  zero  mean;  thus,  we  define  Gq£’  using  (3.19)  for  p  /  0  and  set 
go  =  0.  The  Fourier  coefficients  decay  rapidly  and  we  arrive  at  the  separated 
representation, 


(3.20) 


M 

Gr(x)  =  £  Pn 

m=  1 


E 


r3/2 


pez3,P/o  Lm 


3/2 


e-7T2p2/tm  e27rix-p 


In  (3.20)  let  us  find  M  <  M  such  that 

M  77-3/2 


E  - 

m=M-\- 1 


3/2 


,-7T2/t„ 


<  e. 


We  then  have 


M 


(3.21)  Gger(x)  =  E  Pm  E 


T=1  pGZ3 


t' 


t3/2 


M 


,  — 7T2p 2 /tm  g27rix-p 


E 


t3/2 
+3/T  ' 


m=l 


3.5.  Computing  of  Madelung  Potentials.  Let  us  consider  potential  defined  by 
the  lattice  sum 


(3.22) 


Gm(x)  =  Yj 


(”1) 


ni+n2+ri3 


x  +  n 


1 

bd 


E' 

nez3 


^_\yni+n2+n3 

I  lx  +  nil 


describing  the  crystal  lattice  of  NaCl.  The  function  Gm(x)  is  also  the  Green’s 
function 


V2Gm(x)  =  4t t<5(x), 
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in  B  =  [—1/2, 1/2] 3  with  the  zero  boundary  conditions.  The  second  term  in  (3.22) 
evaluated  at  zero  is  the  Madelung  constant, 


p  =  E 

nez3 


/  (—  l)ral+n2+ra3 


A  careful  computation  of  the  value  of  Madelung  constant  fi  =  —1.74756459...  can 
be  found  in  e.g.  [15],  and  this  potential  and  the  constant  have  been  computed  by 
a  variety  of  methods  [25]  making  it  easy  to  compare  with  our  approach.  We  have 
from  (3.15)  and  (3.16)  an  approximation 

(3.23) 

M  M 

G°Mr  (x)  =  TTUT  -  E  fr>e"iml|X|12  +  E  PmQ(tm,Xi)Q(trn,X2)Q{tm,X3 ), 


where 

Q(t,z)  =  £(- ire-^+")2 

nEZ 


y^(e~4T^2+n)2 

nEZ 


e-4r(x/2+n+l/2)2^ 


For  the  Madelung  constant  we  have 
(3.24) 

M  MM 

E'e“tml|n||2(— 1)"1+"2+"3  =  Eft»E (-ire-*-"2]3  -  Ea 

m— 1  nEZ3  ra— 1  nEZ  m=  1 


Let  us  compute  the  Fourier  coefficients  for  Q(r,  x),  a  periodic  function  with  the 
period  2.  We  have 


qk(T)  =  \  J  Q(r,x)e  mxkdx , 

and,  extending  integration  to  (— oo,+oo), 

<7fc(r)  =  r  (e"4™2  -  e~AT(x+1/2)2)e-2™kdx  =  EL  e-fe2-2/(4r)(1  _ 
J— oo  2-yT 

so  that 


Q(Tai)  =  E  ~r 


/L  „-(fe+l/2)27r2/rp27i-ia;(fc+l/2) 


fee  z 


Thus,  we  arrive  at  (3.23)  and  an  expression  for  the  Madelung  constant, 
m  1  M 

\x  =  7T3/2  E  /WE  e“(n+1/2)2,r2/‘m]3  -  E  pm  =  -1.7475645946... 

m— 1  nEz 


We  note  that  computing  of  slowly  (conditionally)  convergent  sum  has  been  reduced 
to  an  explicit  expression  involving  only  a  small  number  of  fast  convergent  sums  over 
a  one-dimentional  lattice. 

We  note  that 


1 

"w 


M 

Epme-*-||x|1 

m— 1 


<  e 


outside  the  ball  of  radius  1/2  and,  thus,  on  the  boundary  of  the  cube  B.  Also 
Q(i,±l/2)  =  0  for  any  t  >  0  so  that  G^pr(x)  is  the  Green’s  function  of  the 
homogeneous  boundary  value  problem  in  B  with  accuracy  e. 
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Figure  3.2.  A  periodic  part  of  the  external  contribution  to  the 
Madelung  potential  (a  2D  slice  at  z  =  0). 


4.  Results  for  Task  2:  Representation  via  Gaussians  for 

TWO-PARTICLE  PROBLEMS 


One  of  the  key  aspects  in  our  approach  is  to  represent  cusps  due  to  electron- 
electron  interaction.  For  illustration,  we  consider  the  simplest  two-particle  problem, 
the  Helium  atom.  We  have  the  Hamiltonian, 


(4.1) 


2  1 

||x2||  +  ||xx  —  x2  II  ’ 


and  need  to  solve  the  eigenvalue  problem, 


(4.2) 


HV  =  EV. 


Let  us  seek  the  ground  state  wave  function  in  the  form 


(4.3)  $(X1,X2)  =  ^]1/;me-«-llxil|2-/3-llxi-^i|2-Tml|x2||2^ 

In  quantum  chemistry  such  form  of  the  wave  function  has  been  used  before,  see  e.g. 
[37],  where  the  exponents  am,  j3m,  7m  were  selected  upfront  and  the  coefficients 
ipm  computed  to  minimize  the  energy. 

In  our  approach  we  solve  for  the  exponents  am,  (3m,  7 m ,  the  coefficients  ipm 
and  the  number  of  terms  so  that  the  wave  function  satisfies  the  equation  with  a 
prescribed  accuracy.  Such  approach  is  feasible  since  we  can  show  that,  after  we 
rewrite  (4.2)  as  an  integral  equation  (similar  to  what  we  do  in  [28,  30,  50,  51]), 
the  form  of  the  wave  function  (4.3)  is  preserved  under  the  iteration  performed  to 
solve  the  equation.  Since  each  iteration  increases  the  number  of  terms,  the  problem 
then  is  to  develop  an  algorithm  for  reducing  the  number  of  terms  in  the  separated 
representation  (4.3). 
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This  leads  to  the  following  problem:  assuming  that  we  are  given  a  function 

M 

(4.4)  f(x  t,X2,...xN)=  y  wme~ 

m= 1 

where  Xj  €  [0,1],  Tmj  >  0,  and  wm  >  0,  find  a  shorter  optimal  representation  of 
the  same  form.  Namely,  if  the  number  of  terms  in  (4.4)  is  not  optimal  for  a  given 
accuracy  e,  find  g, 

M 

(4.5)  g(x i,x2,...xN)  =  y  wrrae~^=1'rm’jXN 

m= 1 

11/  —  g\\  <  e,  such  that  M  <  M. 

1.  We  have  solved  this  problem  for  N  =  1  in  [14]  and  have  developed  a  pre¬ 
liminary  algorithm  for  N  =  2.  We  now  understand  analytic  and  geometric 
ingredients  of  this  problem.  We  computed  optimal  representations  in  a  num¬ 
ber  of  two  and  three  dimensional  examples  and,  currently,  we  are  developing 
further  additional  components  of  the  algorithm.  For  example,  for  dimen¬ 
sions  greater  than  two,  we  need  to  group  the  initial  exponents  according  to 
their  spatial  distribution,  and  then  solve  several  subproblems.  Our  results 
for  N  =  2  should  lead  to  many  applications  beyond  the  scope  of  this  project 
(e.g.,  generalized  Gaussian  quadratures  for  exponentials  for  2D  domains,  thus 
extending  results  in  [12]). 

2.  We  have  developed  a  new  “suboptimal”  reduction  algorithm  and  computed 
several  examples  for  dimension  N  =  3.  One  of  the  examples  is  the  numerical 
construction  of  an  approximation  of  the  Coulomb  Green’s  function.  This  case 
is  of  interest  for  both  Task  1  and  Task  2.  Such  Green’s  function  (which  is  not 
a  convolution)  has  an  analytic  expression  known  since  1960s  providing  for  us 
a  mechanism  to  verify  our  results.  We  note  that  the  analytic  expression  itself 
(e.g.  [34], [40])  has  not  been  used  in  computations. 

3.  We  have  now  a  mathematical  framework  (with  some  proofs)  for  the  construc¬ 
tion  of  the  multiparticle  Green’s  functions  for  the  so-called  Rollnik  classes  of 
potentials  (see  [5]).  Such  potentials  (or  their  approximations)  are  sufficient 
for  computations  with  bound  states.  This  is  important  on  a  conceptual  and 
on  a  practical  level  (e.g.,  we  obtain  convergent  algorithms). 

4.  Our  approach  is  applicable  to  multi-particle  computations  (in  contrast  to 
using  the  Gross-Pitaevskii  equation)  for  Bose-Einstein  Condensate  (BEC). 

We  now  provide  a  more  detailed  description  of  the  results. 

4.1.  Green’s  functions  for  a  central  potential.  For  the  Hamiltonian 

n  =  -\\72~V{r), 

let  us  consider  the  Green’s  functions 

G(p)  =  (H  -  gl)-1 

G0{g)  =  {-l-S72-gl)-\ 


and 
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Proposition  4.1.  Let  V  be  a  potential  of  the  form 

niiyii)  =  E^e_“fc"y"2’ 

k 

and,  p,  <  infcr(7Y),  where  infa(7i)  is  the  lower  bound  of  the  spectrum  ofTL. 

Then  for  any  e  >  0  and  5  >  0,  <5  <  |  |x  —  y 1 1  <  1,  the  Green’s  function  G(/x,  x,  y) 
has  an  approximation  G(/z,x, y)  ,||G(/z,  x,  y)  —  G(/z,  x,  y)||  <  e/||x  —  y||,  where 


(4.6)  G(/z,x,  y)  = 

m 

Remark.  If  instead  of  the  eigenvalue  problem 

(4.7)  mV  =  AH', 
we  solve 

G(/x)^  =  — tf, 

A  —  /i 

then  (£Jo  —  /x)_1  is  the  largest  eigenvalue,  where  Eq  is  the  eigenvalue  corresponding 
to  the  ground  state  of  (4.7).  Thus  we  have  ||G(/z)||2  =  (Eq  —  m)_1  =  dist(/i,  cr(H)). 
We  can  use  the  power  method  to  find  eigenvalues  and  the  overall  approach  is  simply 
the  inverse  power  method. 


The  Green’s  function  G  satisfies  the  Lipmann-Schwinger  equation, 

G(fj,)  =  Go{fi)  +  G0(fi)VG(n), 

from  which  we  obtain 

g  =  (i-  G0y)_1G0  =  gI/2(i-  gJ/2vgJ/2)_1gJ/2. 

1  /  2  1  /  2 

It  follows  from  Proposition  4.4  in  Appendix  A,  that  the  operator  X  —  G0  VG 0' 
is  positive  definite.  Also,  it  readily  follows  that  the  operator  Gj/2UGj/2  is  positive 
definite  as  well  since  V  is  an  operator  of  multiplication  by  a  positive  function.  This, 
in  turn,  implies  that  || X  —  Gj  2UGj/2||  <  1. 

Following  Proposition  4.2  in  Appendix  A,  we  have 


(4.8) 


G  =  G/2  [|(X+[Go/"PGo/z]"j)Gi 
3=0 


1/2 
0  ' 


Let  us  examine  the  operator  G/2UGj/2.  Using  Proposition  4.5,  we  have 

-i  [•  OO  pOO 

=  ds  i/  ds2  e_3^2(e_2si+e-2S2)+2(Sl+S2) 

J  —  OO  1—00 

"2"2si  P-“fclly!l2  P-lty-z!l2e2S2 


[G2/Vg/2](x,z) 


k 

9-l|x-y|| 


dy- 


Using 

(4.9) 


IV 

/  exph^ajfe -y)2]^ 

j=l 


\iV 
2^7  =  1 


exp[-— ^ ^  a»aj(ari  -  ^j)2], 


l=i 


El\ 

i=i  ai 


i>j 
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and  4.1,  we  obtain 
[G^VG^Kx,  z)  = 


x 


TrS/2 

k 


Vk 


ds  i 


ds- 


e 


(e2 


uk 


i)3/2 


e2sl  +e-*s2  +u 


e2s1  +e^s2  +u 


iP-S-A 

1  e2sl 


E  ufc 

+  e2s2+„ 


This  double  integral  can  be  discretized  using  the  trapezoidal  rule  in  two  dimensions, 
similar  to  the  construction  in  [29,  30,  13],  to  yield  an  approximate  representation  of 
the  form  stated  in  (4.14).  At  this  point  we  are  not  counting  the  necessary  number 
of  terms,  just  noting  that  they  are  of  the  desired  form.  The  powers  of  T0 ,  , 

as  needed  in  (4.8),  will  have  a  similar  integral  representation  with  the  number  of 
integrals  doubling  as  we  step  from  j  to  j  +  1.  Again,  using  the  trapezoidal  rule,  we 
obtain  terms  of  the  desired  form  and,  as  the  number  of  retained  terms  in  (4.8)  is 
finite  for  a  given  accuracy  e,  the  process  will  terminate  after  a  finite  number  of  terms. 
Since  we  use  potentials  without  singularities,  the  norm  G(/x,  x,  y)  —  G(/i,  x,  y)  will 
be  controlled  by  selecting  the  parameters  of  the  trapezoidal  rule.  This  argument 
provides  the  form  of  the  Green’s  function  but  not  the  number  of  terms. 

In  order  to  make  this  computation  practical,  we  turn  to  the  second  part  of 
Proposition  4.2.  We  have  a  recursion, 

yn+i  =  2 yn  -yn(z-  g\/2vg\/2)  yn,  with  =  i, 

where  yn  —>  (1  —  G^VG^2)-1. 

1  /2 

Let  us  use  a  representation  of  the  kernel  G0  via  Gaussians  as  it  follows  from 
Appendix  B, 

gJ'V* -y>  =  £»«-' 

I 


1/2  1/2 

Computing  G0  PG0  ,  we  have  an  approximation 


(4.10) 


E 

LkM 


giVkgv 


0-«il|x-y|! 


a-“fcl|y||2  „-si'l!y-zll 


dy 


Using  (4.9),  we  arrive  at 
(4.11) 


l,k,l’ 


givkgi'{ 


Sl  +  Si>  +  uk ' 


-)3/2 e 


sl+sl,+uk  I 


»'‘l' 

Sl+S,/+Ufc  I 


S(+S,/+Ufe 


INI! 


which  is  of  the  form  in  (4.14).  We  now  need  to  reduce  the  number  of  terms  in  (4.11) 
and,  for  this  reason,  arrive  at  the  following  problem:  given  a  function 

f(x)  =  Y^wke-a*Xl-a*x*-a*X3 , 

k 


approximate  /  by  a  function  /  of  the  same  form  but  with  a  fewer  number  of  terms. 

We  note  that  the  form  in  (4.14)  is  preserved  under  further  iteration,  namely,  if 
we  compute  the  integral  (for  a  single  term), 


(4.12) 


Ilkm  —  giVkWrn 


=  -«l||x-y||  p— “felly  1 1  p-Tm  |  |y  1 1  -CTm||y-z||  -€m|f 


dy. 
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and  use  (4.9),  we  obtain 

(4.13)  Iikm  =  giVkWmi - )3/  e  clkm  Hkm"  (5m  °lfcm 


Qfcra 


where  =  s;  +  Uk  +  am  +  rm.  This  gives  us 


_  Sl(uk+Tm)  ,  _  SjCJm  ~  _  ,  (ttfc  + 

7~lkm  —  5  ®lkm  —  ?  <,lkm  —  swi  '  ; 

Qfem  Qfem  Qfcm 


and 


Wlkm  =  9lVkWm[ 


C-lkm 


-)3/2, 


as  new  parameters. 

4.2.  Green’s  functions  for  multiparticle  systems. 

4.2.1.  Confining  potential.  Consider  the  Hamiltonian  of  the  system  with  K  nuclei, 

JVe  ..  K 

*-»-  &-E  A). 


J  =  1 


fc  =  l 


The  Green’s  function. 


can  be  approximated  by 


G(p)  =  (H-pI)-\ 


(4.14) 

G(/i,xi,yi,x2,y2,...)  = 


for  p  <  Eq,  where  E0  is  the  smallest  eigenvalue  of  H  . 

4.2.2.  H  Hamiltonian  with  the  electron- electron  interaction.  Consider  the  Hamil¬ 
tonian  of  the  system  with  if  nuclei  and  Ne  electrons, 


Ne 


«  =  ir)  +  E 


a: 


1 


i=i 


fc=i 


x7-  -  rfc 


*>J 


The  Green’s  function, 


can  be  approximated  by 


G(p)  =  (H-pI)~\ 


G(/i,  Xi,  yi,  x2,  y2, . . . )  =  ^ 


-  EdEfc  am.fcllxj-rfclP+Efc  bm,fc||yi-rfc||2] 


10™  e  ^LZ^fc 


x  e 


Cm. if  Xi 


y j  1 1  Yli>j  [9m, ij  |  |xi  Xjll  fm,ij  |  |yi  Yj  |  |  ] 


(4.15) 

for  p  <  E0,  where  E0  is  the  smallest  eigenvalue  of  TL  .  We  conjecture  that  it  may 
be  possible  to  simplify, 

G(p,x  i,yi,x2,y2, . . . )  =  ^wm  «m.fc(l|x,--rfc||  +|  — rfc  1 1  ) 


x  e“  cm.wllxi-y^ll  -Z)4>J  Sm,«(llx»-xj||  +||y.-ydl  ) 


(4.16) 
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4.2.3.  Spectral  Projectors.  Representation  of  a  spectral  projector  computed  for  the 
Green’s  function  has  the  structure  as  the  Green’s  function  itself  (unlike  the  eigen¬ 
functions,  which  need  additional  terms  or  factors  to  account  for  momentum  indices). 
We  construct  the  spectral  projector  from  the  sign  function  using 

(4.17)  P^y)  =  Y  ^(x)Vh(y)  =  (X-  sign (H  -  pJ))/2, 

A  j<fi 


where  ifjj  is  an  orthonormal  basis  of  eigenfunctions  and  the  sign  function  is  defined 
on  (—oo,  oo )  by 


(4.18) 


sign(A) 


1  A  >  0 

0  A  =  0 

-1  A  <  0 


If  an  operator  or  its  projection  on  the  discrete  spectrum  are  written  as 

(4.19)  H(x,  y)  =  Y  V-'.bxi'.j(y! 

3 


with  A j  real,  then 


(4.20)  sign(W)(x,y)  =  ^sign(AJ-)V’j(x)^(y). 

i 


4.2.4.  Recursive  construction.  We  use  a  polynomial  recursion  (see  e.g.  [4,  33,  7]) 
to  compute  sign(7d): 

,  .  y0  =  n/\\n\\2 

(4.21) 

34+1  =  (334  -  343)A  k  =  0,1,... 

Other  polynomials  may  be  used  in  place  of  the  one  above;  see  [33]  for  a  discussion 
of  the  various  choices.  It  is  easy  to  demonstrate  that  34-  — >  sign(W)  in  (4.21)  (see 
e.g.  [7]).  The  number  of  iterations  needed  for  (4.21)  to  converge  to  accuracy  e  is 
0(clog2  k  +  log2  log2(l/e)),  where  k  is  the  condition  number  of  34  • 

We  note  that  computation  of  the  projections  via  (4.21)  has  qualitatively  different 
properties  than  that  of  the  direct  computation  of  individual  eigenfunctions.  As  in 
the  case  of  the  Green’s  function  the  iteration  in  (4.21)  does  not  change  the  form  of 
the  Gaussian  representation. 


4.2.5.  Bound  state  solutions  of  the  Lippman- Schwinger  Equation.  Consider  a  self- 
adjoint  operator 

(4.22)  H  =  H0-V 

and  the  Green’s  function 

g(z)  =  {H-zi)-\ 

for  z  cr(H),  where  a(Tt)  denotes  the  spectrum  of  the  operator  H.  Assuming 
2  <r(7io)  and  introducing 

Go(z)  =  (Ho-zI)-1, 

we  (formally)  have  the  Lippman-Schwinger  equation  for  Q, 

Q  =  Go  +  GoVG, 

or  G  =  (1  —  GoV)~1Go-  Alternatively,  we  can  write 

G  =  Gl/2(X  -  Gl/iVGl,2)-1Gl,i . 
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To  construct  (I  —  Gq^VQ^2)  1,  we  will  use 

Proposition  4.2.  1.  IfB  is  a  bounded  positive  definite  operator,  then  its  inverse 

has  a  product  representation 

OO 

(4.23)  B~x  =  -  nBf], 

j= o 

where  k  <  1/||£>||. 

2.  The  iteration 

yn+i  =  2yn  -  ynByn,  with  y0  =  «x, 

converges  quadratically  to  the  inverse,  yn  — >  B -1  as  n  — >  oo  and 

n 

yn+i  =  k  n  [J  +  (J  -  kS)2'L  n  =  o,  1, . . . . 

j= o 

Proof.  We  have  .EU1  =  k(I  —  (X  —  k13))_1.  Since  B  is  a  bounded  positive  definite 
operator,  selecting  k  <  1/||£>||  implies  that  \\T  —  kB\\  <  1.  Using  the  product  form 
of  the  converging  series  B^1  =  n[l  +  (I  —  kB)  +  (I  —  kB)2  +  ...],  we  obtain  (4.23). 
The  iteration  simply  generates  the  same  product  as  in  (4.23).  □ 

Proposition  4.3. 

1.  Let  B  be  a  bounded  operator.  If  the  inverse  operator  £>_1  exists,  then  it  has 
a  product  representation 

OO 

(4.24)  B~1=nY[[l+{l-KB*B)2i]B*, 

3=0 

where  k  <  1/||£>*  B\\. 

(a)  The  iteration 

yn+ 1  =  2yn  -ynByn,  with  ;y0  =  kB\ 

converges  quadratically  to  the  inverse,  yn  — >  B as  n  — >  oo  and 

n 

yn+ 1  =  k\\{1+(1  -  nB*B)2j]B*,  n  =  0,1,.... 

3=0 

Proof.  We  have  B~x  =  (B*B)~1B* .  Since  B*B  is  a  bounded  positive  definite  oper¬ 
ator,  we  use  Proposition  4.2  to  compute  (B*B)~1.  □ 

Proposition  4.4.  If  the  potential  V  is  in  the  Rollnik  class  and  z  (7(74),  then 
B  =  T  —  G\J2  (z)VGq^2  (z)  is  a  bounded  operator.  If  z  el  and  z  <  min(cr(74)),  then 
B  is  a  positive  definite  operator  and  \\Gq^  {z)VGq^  (z)\\  <  1- 

Proof.  Let  us  consider  the  inner  product  (Bx,  x)  =  (x,  x)  +  (VGq/'2x,  Qq2x').  Setting 
y  =  Q\I2x,  we  have 

(Bx,  x)  =  ( Go1/2y ,  Gq1/2V )  +  ( Vy ,  y)  =  (Gffly,  y)  +  (' Vy ,  y)  =  ((74  -  zl)y,  y)  >  0 

for  any  x.  Since  G^2  (z)VG^2  (z)  is  also  positive  definite  implies  that  ||  (^)V^/q//2  (^)  || 

1.  □ 
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Using  Proposition  we  have  that  if  B  is  bounded,  then 


(4.25)  G  =  kQI/2  \{(X  +  [(1  -  k)1  +  kG^VGI'Y  )  Gl/2. 

j= o 

If  \\G^'2 (z)VGY (z)\\  <  1  then  ll#l!  <  1  and  we  choose  k  =  1  to  obtain  from  (4.25) 

OO 

(4.26)  G  =  Go/2  I](J  +  [GTvGTfYJ2, 

3=0 

an  alternative  of  the  usual  Born  series. 

4.2.6.  Integral  Representation  and  Approximation  of  G o  •  In  problems  of  quantum 
mechanics  we  choose 

7Yo  =  4v2, 

and  construct  an  approximation  of  the  kernel  of  Qq,  a  >  0,  where 

0o(A)  =  (Ho  —  AT)-1. 

We  first  recall  the  kernel  of  etv2,  the  heat  kernel, 

(4.27)  A'etV2(x- y)  =  (47rt)_3e_3llx~yll2t  \ 

and  the  integral  representation  of  the  positive  powers  of  a  self-adjoint  positive 
definite  operator, 


L~a  = 


1 


P(a) 


e~tLta~1dt. 


We  have  for  A  =  —  ^  \r  <0  , 

1 


^  oq  =  (-aV2  +  -mV)-“  = 


T(a) 


etv  e_M  ttOL~1dt. 


r(a)(47r)2  Jo 


o-^te-Mx-yf^-ldt. 


Using  (4.27),  we  have 

(4.28)  £o(x-y)  = 

Since 

r°°  212  h  h 

(4.29)  /  e~a  t^dt  =  21"C(-)CK c(ab)  =  21“C(-)CK _c(o6), 

Jo  a  a 

where  Ka  is  the  modified  Bessel  function  of  the  second  kind.  This  formula  follows 
from  [27,  Formula  8.432(6)]  , 

roo 

K c(z)  =  2c~1z~c  /  e-'e-®2  U'df,  for  |  arg^l  <  3?z2  >  0. 

Jo  2 

Thus,  Go  is  a  radial  function  that  can  be  explicitely  described  as 


r(a)7T2  r ' 


where  r  =  |x|.  In  particular, 


Go  (x) 


I  e— mIx( 
277  Ixl 
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and 


To  approximate  Qq  as  a  sum  of  Gaussians,  we  change  variables  t  =  e  2s /4  in  (4.28) 
to  obtain 

Proposition  4.5.  The  kernel  of  Qq  has  an  integral  representation 


ol-a  poo 

(4.30)  So(x-y)  =  zrrw— 3-  / 

1  (a)7r2  J-oo 


e-\\K-yW*e*‘ e-\^e-*°+(Z-2a)sds 


For  a  =  1  we  obtain  from  (f-30)  a  familiar  representation  (see  [29,  30,  13],), 
Go(x-y)  =  -4  I"  e-^-y^se-^2e~2a+sds, 

7T  2  J  —  oo 

and,  for  a  =  1/2,  an  integral  representation  for  the  kernel  of  , 


£o1/2(x-y)  = 


23 


-||x-y||2e2se-Ip2e-2s+2S(is^ 


5.  Results  for  Task  3:  Multiparticle  Schrodinger  Equation 


This  part  of  the  project  is  the  beginning  of  the  program  to  develop  accurate 
methods  for  solving  equations  of  multiparticle  quantum  mechanics.  In  this  section 
we  will  descibe  what  was  accomplished,  compare  with  existing  state  of  the  art,  and 
indicate  directions  of  further  development. 

Our  goal  was  specifically  to  address  the  issue  of  size-consistency.  We  first  devel¬ 
oped  a  “center-of-mass”  principle  upon  which  an  algorithmic  size-consistent  method 
can  be  built.  We  were  then  advised  that  the  inter-electron  cusp  was  a  crucial  is¬ 
sue,  and  so  we  developed  a  principle  upon  which  a  method  capturing  this  cusp  can 
be  built.  As  these  methods  were  becoming  large  and  complex,  we  went  back  and 
worked  on  the  details  of  the  basic  method  upon  which  they  were  grown,  namely  the 
use  of  separated  representations  for  wavefunctions.  In  the  next  section  we  give  a 
technical  introduction  to  the  multiparticle  Schrodinger,  and  the  outline  of  how  sep¬ 
arated  representations  can  be  used  to  approximate  the  wavefunction.  The  following 
section  sketches  the  method  for  size-consistency,  and  the  final  section  comments  on 
the  inter-electron  cusp. 


5.1.  Technical  introduction;  Approximating  the  Wavefunction  with  an 
unconstrained  sum  of  Slater  Determinants.  The  multiparticle  Schrodinger 
equation  is  the  basic  governing  equation  in  quantum  mechanics.  We  consider  the 
time-independent  case,  and  fix  the  nuclei  according  to  the  Born-Oppenheimer  ap¬ 
proximation,  so  the  equation  describes  the  steady  state  of  an  interacting  system  of 
electrons.  For  each  of  the  N  electrons  in  the  system  there  are  three  spatial  vari¬ 
ables  r  =  ( x,y,z ),  and  a  discrete  spin  variable  er  taking  the  values  {—  |},  which 

we  combine  and  denote  (r,  a)  by  7.  The  Hamiltonian  operator  74  is  a  sum  of  a 
kinetic  energy  operator  T,  a  nuclear  potential  operator  V,  and  an  electron-electron 
interaction  operator  W,  defined  by 


N 


N 


(5.1) 


n  =  r  +  v  +  w  = 


N  N 


*=1 


1 

r,;  -  r., 
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where  A  $  is  the  three-dimensional  Laplacian  acting  in  the  variable  r,  and  v(r)  is  a 
sum  of  terms  of  the  form  za/ ||r  —  Ra||  from  a  nucleus  at  position  Ra  with  charge 
za-  The  antisymmetric  eigenfunctions  of  H.  represent  electronic  states  of  the  system 
and  are  called  wavefunctions.  Antisymmetric  means  that  under  the  exchange  of  any 
two  coordinates,  the  wavefunction  is  odd,  e.g.  ip( 72, 7i,  ■  •  ■ )  =  —'0(71,72,  •  •  •  )•  The 
bound-state  wavefunctions  have  negative  eigenvalues,  and  are  of  greatest  interest,  so 
we  will  focus  on  the  wavefunction  with  the  most  negative  eigenvalue.  In  summary, 
our  goal  is  to  find  the  most  negative  (discrete)  eigenvalue 

(5.2)  Hip  =  A  ip , 


subject  to  the  antisymmetry  condition  on  ip. 

Analytic  methods  can  give  qualitative  results  about  its  solutions,  and  determine 
limiting  cases,  but  most  quantitative  results  must  be  obtained  numerically.  Al¬ 
though  the  equation  is  a  ‘simple’  eigenvalue  problem,  its  numerical  solution  presents 
several  serious  difficulties,  among  them  the  large  number  of  variables  and  the  an¬ 
tisymmetry  condition  on  the  solution.  The  simplest  method  that  addresses  these 
two  difficulties  is  Hartree-Fock  (HF),  which  uses  the  antisynnnetrization  of  a  single 
product  to  approximate  the  A-particle  wavefunction,  i.e. , 

N 

(5.3)  V’HF  =  Aj^^i(7i) . 

i= 1 


Any  approximation  ^  to  the  wavefunction  ^  can  be  substituted  into 


(5.4) 


W,0) 

M 


to  obtain  an  upper  bound  on  the  lowest  value  of  A  that  solves  (5.2).  Substituting 

(5.3)  into  (5.4),  one  can  derive  a  system  of  equations  for  (pi  to  minimize  (5.4).  The 
resulting  0hf  will  best  approximate  ip,  in  the  sense  of  providing  the  best  estimate 

(5.4) .  _ 

To  improve  upon  HF,  it  is  natural  to  consider  a  sum  of  products 

r  N 

(5.5)  V’(r)  =  Py'(7i)  • 

1=1  i—1 


The  coefficients  Si  are  not  strictly  necessary,  but  they  allow  us  to  assume  ||0(||  =  1. 
Many  methods  are  based  on  this  form,  and  the  distinction  is  in  how  they  use  it.  The 
Configuration  Interaction  (Cl)  method  (see  e.g.  [48])  chooses  the  functions  (p\  from 
a  preselected  master  set  of  orthogonal  functions  and  decides  on  a  large  number  r 
of  combinations  to  consider,  based  on  excitation  level.  Substituting  (5.5)  into  (5.4) 
leads  to  a  matrix  eigenvalue  problem  that  can  be  solved  for  the  scalar  coefficients 
si .  The  Multi-Configuration  Self-Consistent  Field  (MCSCF)  method  (e.g.  [24,  17]), 
chooses  a  pattern  of  excitations  similar  to  Cl,  but  then  solves  for  the  master  set 
of  orthogonal  functions  as  well  as  the  scalar  coefficients.  Many  variations  and 
combinations  of  these  methods  have  been  developed,  and  indeed  there  is  a  whole 
industry  in  producing  them. 

We  demonstrate  a  method  that  also  uses  a  wavefunction  of  the  form  (5.5)  but 
without  constraints  such  as  orthogonality  on  the  </>(.  By  removing  these  constraints 
we  produce  much  better  approximations  at  much  smaller  r  than  existing  methods 
allow.  In  another  context  [11]  we  have  given  examples  where  removing  constraints 
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produces  expansions  that  are  exponentially  more  efficient,  i.e.  r  =  N  instead  of  2N 
or  r  =  log  N  instead  of  N .  For  example,  in  our  approach  we  can  have  a  two-term 
representation 

0i (71)^2(72)  •  •  •  0jv(7iv) 

(5.6)  +c  [0i (71)  +  0iv+i (71)] [^2(72)  +  ^+2(72)]  •  •  •  [</>jv(7at)  +  <!>2n{in)\  , 

where  form  an  orthonormal  set.  To  represent  (5.6)  while  requiring  all 

factors  to  come  from  a  master  orthogonal  set  would  force  one  to  multiply  out  the 
second  term  and  thus  obtain  a  representation  with  2N  terms.  It  is  common  sense 
that  removal  of  constraints  could  produce  better  results,  and  steps  in  that  direction 
have  been  taken  (e.g.[46,  1,  23,  19,  20,  52,  42]).  These  works,  however,  were  only 
able  to  partially  remove  the  constraints,  and  so,  we  claim,  did  not  achieve  the  full 
potential. 

We  will  use  a  Green’s  function  iteration  to  move  a  trial  wavefunction  toward  the 
minimum  of  (5.4)  without  using  (5.4)  directly.  This  iteration  was  introduced  in 
[36,  35]  and  used  in  e.g.[30].  Define  the  Green’s  function 

(5.7)  =  (T  -  pi)-1 , 
for  p  <  0.  The  Green’s  function  iteration  is 

9n  =  -S„n[(V  +  W)/n] 

(5.8)  pn+ 1  =  pn-((V  +  W)fnJn~gn)/\\gn\\2 

fn+ 1  —  9n/\\9n  ||  • 

The  Green’s  function  iteration  is  essentially  an  inverse  power  method.  The  conver¬ 
gence  rate  is  only  linear,  but  if  the  initial  p  can  be  chosen  near  to  but  less  than 
the  lowest  eigenvalue,  then  the  error  will  decrease  by  a  substantial  fraction  at  each 
iteration,  and  not  many  iterations  will  be  needed.  We  use  /  to  denote  the  number 
of  Green’s  function  iterations  needed. 

We  use  approximate  wavefunctions  of  the  form  (5.5),  with  r  fixed.  The  iteration 

(5.8)  does  not  directly  produce  an  approximation  of  the  same  form,  so  we  modify 
it  by  defining  gn  to  be  the  function  of  the  form  (5.5)  that  minimizes 

(5-9)  \\9n-(-G^n[(V  +  W)fn])\\. 

In  order  to  assure  convergence  to  an  antisymmetric  solution,  we  use  the  pseudo¬ 
norm  induced  by  the  pseudo  inner  product  (•,  -)^t  =  (A(-), -4(-)),  as  we  did  in  [11]. 
Constructing  gn  is  the  most  challenging  part  of  the  method,  and  requires  the  bulk 
of  our  effort.  To  simplfy  notation,  we  now  suppress  the  iteration  index  n  and  set 
0  =  fn  and  0  =  gn- 

We  begin  with  some  approximation  0  (such  as  0  itself)  and  will  iteratively 
improve  it.  The  outermost  loop  of  our  iteration  is  simply  to  repeat  our  refinement 
until  it  appears  that  0  has  converged.  For  the  computational  cost  estimates  we 
denote  the  number  of  repetitions  by  K .  To  refine  our  representation  we  loop  through 
the  variables  (electrons) .  The  functions  in  variables  other  than  the  current  variable 
are  fixed,  and  the  functions  in  the  current  variable  will  be  modified  to  minimize  the 
overall  error  ||0  —  0||a-  This  Alternating  Least-Squares  (ALS)  approach  is  well- 
known  (see  e.g.  [31,  39,  41,  16,  18,  49]).  We  will  alternate  through  the  directions, 
but  for  ease  of  exposition  we  describe  the  k  =  1  case.  So,  (j>lk  is  fixed  for  k  >  1,  and 
we  will  solve  for  the  values  of  for  all  l. 
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To  refine  in  the  current  variable,  we  set  up  and  solve  a  linear  least-squares 
problem.  The  normal  equations  for  a  least-squares  problem  are  derived  by  taking  a 
gradient  with  respect  to  the  free  parameters  and  setting  this  equal  to  zero.  As  long 
as  if)  is  linear  and  not  degenerate  in  these  parameters,  the  resulting  equations  are 
linear  and  have  a  unique  solution.  Usually  these  free  parameters  are  coefficients  of 
the  representation  in  some  basis.  We  instead  take  the  parameters  to  be  the  point 
values  of  our  functions  (f>[ ,  or,  formally,  as  the  coefficients  of  the  point  evaluation 
functional  (7) .  The  formulas  that  we  derive  can  be  used  with  a  fixed  basis,  but  are 
stated  independent  of  the  basis.  We  still  obtain  linear  normal  equations 

(5.10)  Ax  =  b  , 

but  now  b(l)  is  a  function  of  7,  x(l')  is  a  function  of  7',  and  A(l,l')  is  an  integral 
operator  mapping  functions  of  7'  to  functions  of  7.  The  kernel  of  A  is  defined  by 

IN  N 

(5-ii)  A(i,i')(  7,7')  = 

\  i= 2  i—2 

where  the  point  evaluation  functionals  are  acting  in  the  i  =  1  direction.  The 
functions  in  b  are  defined  by 

r  /  N  N 

(5.12)  6(0(7)  =  E  +  W]  n  n  & 

m  \  i—1  i—2 

Once  A  and  b  have  been  constructed,  we  will  apply  the  Conjugate  Gradient  iterative 
method  (see  e.g.  [26])  to  solve  (5.10).  We  initialize  with  r  =  b  —  Ax,  v  =  r,  and 
c  =  (r,  r),  and  then  the  core  of  the  method  is  the  sequence  of  assignments  z  <—  Av, 
t  <—  c/(v,  z),  x  <—  x  +  fv,  r  <—  r  —  tz,  d  <—  (r,  r),  v  <—  r  +  (d/c)v,  and  c  d,  applied 
iteratively.  We  use  S  to  denote  the  number  of  conjugate  gradiant  iterations  needed. 
Thus  x  is  constructed  using  only  matrix-vector  products  and  vector  additions,  all 
which  are  compatible  with  our  formulation  with  integral  operators.  The  conjugate 
gradient  method  applies  only  to  positive-definite  operators.  Our  operator  A  is  only 
semidefinite  due  to  the  nullspace  in  the  antisymmetric  pseudonorm.  Fortunately,  b 
was  computed  with  the  same  pseudonorm  and  has  no  component  in  the  nullspace 
of  A. 

One  advantage  of  using  this  iterative  method  with  integral  operators  is  that  our 
algorithm  is  “basis-free”.  The  representation  of  x  can  naturally  be  adaptive  in  7, 
for  example  refining  near  the  nuclei  as  indicated  by  the  refinement  in  b.  For  the 
estimates  of  computational  cost,  we  use  M  to  denote  the  cost  to  represent  a  function 
of  7,  or  integrate  such  a  function.  The  antisymmetry  contraint  requires  N  <  M, 
and  in  general  we  expect  M  to  be  much  larger  than  N.  For  our  numerical  results, 
we  use  adaptive  polynomial  multiwavelets,  following  [28,  30].  In  those  works  it  was 
shown  that  this  basis  effectively  eliminates  basis-set  error  within  HF. 

We  are  left  with  the  problem  of  how  to  construct  A  in  (5.11)  and  b  in  (5.12).  We 
haved  develop  the  machinery  and  algorithms  for  computing  these  antisymmetric 
inner  products.  Our  formulation  uses  low-rank  perturbations  of  matrices,  thus 
avoiding  cofactor  expansions. 

The  computational  cost  for  the  whole  method  is  acceptable.  As  noted  above,  the 
cost  depends  on  IV,  r,  M,  /,  K,  and  S.  Although  S  in  theory  could  be  as  many  as 
rM,  we  have  a  very  good  starting  point,  and  so  expect  only  a  very  small  constant 
number  to  be  needed.  We  use  M\ogM  to  denote  the  cost  to  convolve  a  function 
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of  7  with  1/ 1| r || .  Some  Poisson  solvers  achieve  this  complexity,  but  this  cost  may 
vary  with  the  choice  of  basis.  We  use  L  to  denote  the  number  of  terms  used  to 
approximate  the  Green’s  function  with  Gaussians.  The  final  computational  cost  is 
then 

(5.13)  0(KIr2N2[L(N  +  MlogM)  +  S(N  +  M)]). 

For  comparison,  the  cost  to  evaluate  a  single  instance  of  Lowdin’s  rules  is  0(N2(N+ 
M)).  The  size  r  needed  in  practice,  and  how  it  depends  on  the  various  parameters 
in  the  problem,  is  still  an  open  question. 

We  are  working  on  initial  implementation  to  verify  the  correctness  of  the  formulas 
and  algorithms  that  we  have  developed. 

5.2.  Algorithmic  Size-Consistency.  Consider  the  situation  where  our  system 
consists  of  K  non-interacting  (i.e.  well-separated)  subsystems,  each  with  N  elec¬ 
trons.  Suppose  that  each  subsystem  has  a  separated  representation  ipi.  The  wave- 
function  for  the  entire  system  is 

K 

ip « n^- 

j= i 

Its  inherent  complexity  grows  linearly  with  K.  If  each  ipi  has  separation  rank  r,  then 
in  order  to  represent  the  overall  wavefunction  in  the  separated  representation  we 
have  to  multiply  out,  and  so  obtain  rK  terms.  Thus,  this  representation  is  not  size- 
consistent.  (The  terms  “size-extensivity”  and  “size-consistency”  are  often  confused 
in  the  literature,  but  it  appears  that  we  actually  address  “size-consistency”.)  For 
a  fixed  accuracy,  we  want  the  computational  complexity  to  scale  (nearly)  linearly 
with  system  size.  In  classical  particle  systems,  such  scaling  can  be  achieved  via 
organizing  particles  into  hierarchical  groups  and  computing  the  interaction  between 
these  groups  via  e.g.,  multipole  expansions.  In  the  simplest  case,  this  amounts  to 
replacing 


Active  Distant 


Summary 

The  notion  of  groups  is  fairly  straightforward  in  quantum  mechanics,  and  in  the 
non-interacting  case  we  can  simply  not  multiply  out  the  component  wavefunctions, 
and  thus  obtain  a  size-consistent  representation.  Non-interacting  systems  are  not 
very  interesting  of  course,  but  they  give  us  the  inspiration  to  use  the  form 

;=i j= l 

for  interacting  systems.  One  could  plug  this  back  into  itself  and  do  more  levels. 

In  quantum-mechanical  systems,  the  interaction  between  particles  is  much  more 
complicated,  mainly  due  to  the  antisymmetry  constraint  which  allows  “exchange” 
of  particles.  The  interactions  are  thus  nonlocal  and  nonlinear,  so  a  straightforward 
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generalization  of  the  approach  to  classical  particle  systems  is  not  available.  In 
particular,  if  we  have  a  geometry  like 

JO? 

Active  Neighbor  Distant 

then  we  cannot  compute  the  effect  of  the  distant  group  on  the  active  group  without 
accounting  for  exchange  via  the  neighbor  group.  Within  our  ALS  method,  the 
key  ingredient  for  achieving  size-consistency  is  to  have  the  marginal  cost  of  each 
antisymmetric  inner  product  be  independent  of  the  total  number  of  groups,  so  that 
the  total  cost  scales  (nearly)  linearly  in  K.  This  goal  can  be  accomplished  if  the 
effect  of  distant  groups  can  be  summarized  and  then  re-used. 

We  now  demonstrate  how  to  summarize  in  the  case  of  three  groups,  by  computing 


(£  si1 )  (£  ^ )  (£  *3  ).(£  *1 )  (£ ^ )  (£  *!? ) 

l 1  l“2  I3  ^  ^  ^3 


A 


We  use  d*)1  to  denote  the  l\  term  in  group  1,  and  so  on.  We  define  LfTj1 ,  f1)1 )  to  be 
the  matrix  of  inner  products  of  the  single  electron  functions,  as  used  in  Lowdin’s 
rules.  Using  Lowdin’s  rules  and  our  assumption  on  supports,  we  have 


£££ 


L(^iS<1>22)  0 

L(^V^22)  L(^22^33) 

0  L(l>f,d>'3) 


The  off-diagonal  blocks  should  have  rank  approximately  the  number  of  chemical 
bonds  between  the  groups.  To  demonstrate  the  principle,  we  will  suppose  that  the 
rank  is  one.  Applying  the  inverse  of  the  diagonal  blocks,  we  obtain 


£££  |L(l£<^)||L($: 


u. 


,®33)l 

I 

Jil'ih  ( 
u12  \v12  ) 

0 

111 2^2  ( -.J1 

■21  Vv21  ) 

l 

/  hl2h  \* 

u23  \v23  ) 

0 

I2I3I3  (  I2I3I3 \* 
u32  Vv32  ) 

1 

We  now  transform  the  determinant  using  the  following  lemma. 

Lemma  5.1.  (Determinant  of  a  Low-Rank  Perturbation  of  the  Identity)  Let  {ug}^=1 
and  {vglljL-L  be  two  sets  of  vectors.  Then 


Q 

1  +  £  U9vg 


9=1 


l+v*ui  vCu2 

v^Ui  1  +  v5u2 

VQU1  VQU  2 


V1UQ 

V2UQ 

1  +  vQuQ 
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Applying  this  lemma  the  final  determinant  becomes 


1 

/  Z1Z2Z2  \  5tZiZp: 
V  V21  )  U12 

2  0 

0 

(  hlih\*  hW2 

V  v12  )  U21 

1 

0 

\v32  )  U21 

V  v12  )  u23 

0 

1 

/  Z2Z3Z3W  Z2Z2Z3 
\v32  )  u23 

0 

0 

/  Z2Z2Z3  \  *  Z2Z3Z3 
Vv23  )  u32 

1 

which 

can  be  expanded 

as 

1  (v 

Z1Z2Z2  \  *  Z1ZP2 

21  )  U12 

1 

/  ^2^3  Z3  \  *  Z2Z2Z3 

\v32  /  u23 

0 

ZiZ^x*  ZiZ2Z2 
^12  )  U21 

1 

/,  _hl2h\*  Z2Z3Z3 
\v23  )  u32 

1 

//  ZiM2V 
~((V12  ) 


Z  2  Z  2  Z  J 
23 


U. 


)((’ 


hhl'o ' 


U 


Z  i  Z  i  Z ; 


,)((V23i,S 


u^3,»)((v. 


Z2Z3Z3  \  *  Z1Z2Z2  \ 
32  )  U21 


Inserting  this  expansion  and  rearranging,  the  first  term  yields 


^2  /FvZ2 


Zi,Z 


1»*1 


/  ZiZiZ2x*„ZiZ2Zi 
\v12  /  U21 


/__ZiZ2Z2\*  __ZiZ1Z2 
Vv21  1  u12 


and  the  second  term  yields 


-EE  ((v*^)*^2 )  ^ )  1  )  1 

Zi  ,Z^  Z2  5Z2 


^2 Z2 Z3  \  *  _Z2Z3Z3 
vv23  /  u32 


/  Z2Z3Z3  w  Z2Z2Z3 
Vv32  /  u23 


^^u23^3  ((^23^^3  )*^323^  )  |lb(^3^ ,  )  |  (v!^3^)*  |  uf1^2 

\h,l'3 


21 


In  both  cases  the  sum  over  group  three  can  be  performed,  resulting  in  summary 
quantities  that  play  the  role  of  multipole  expansions.  Schematically  we  have 

J 

Active  Neighbor  Summaries 

Given  an  “active”  group,  these  averaged  quantities  and  quantities  associated  with 
its  immediate  neighbors  are  all  that  is  needed  to  compute  the  wavefunction  for 
that  active  group.  The  derivation  of  these  quantities  and  their  use  is  complete  (for 
simple  cases),  and  we  have  begun  the  implementation  of  relevant  algorithms. 


5.3.  Interelectron-cusps.  We  note  that  the  representation  (5.5)  does  not  account 
for  the  inter-electron  cusp  (see  e.g.[47,  43,  38,  44,  45,  37]).  As  with  Cl  methods,  we 
may  still  be  able  to  achieve  small  error  in  the  energy  difference  of  two  configurations, 
which  is  often  the  quantity  of  interest  in  Chemistry.  To  achieve  high  accurary  in  the 
wavefunction  we  are  developing  an  extension  to  (5.5)  that  incorporates  the  cusp. 
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We  consider  wavefunctions  represented  in  the  form 

p= o 

and  have  developed  many  of  the  algorithms  needed  to  use  them.  The  algorithms 
needed  are  quite  complex,  however,  so  we  are  deferring  their  development. 

The  size-consistent  algorithms  will  eventually  need  to  incorporate  the  effect  of 
cusps  as  well,  so  that  we  can  achieve  high-accuracy  and  size-consistency  at  the  same 
time. 


N 


5Zwp(Hri-rj||) 


e  n 


J= 1 k- 1 


6.  Papers  in  preparation 

As  we  have  obtained  a  large  number  of  new  results,  we  have  several  papers  in 
preparation.  We  list  them  here. 

•  Tentative  title:  “Multiresolution  separated  representations  of  lattice  sums”, 

G.  Beylkin,  L.  Monzon  and  R.  Harrison. 

(Multiresolution  separated  representations  of  Green’s  functions  satisfying  bound¬ 
ary  conditions  and  associated  fast  algorithms). 

•  Tentative  title:  “Separated  representations  of  lattice  sums  for  oscillatory  Green’s 
functions”,  G.  Beylkin,  L.  Monzon  and  R.  Harrison. 

(Separated  representations  for  oscillatory  Green’s  functions). 

•  Tentative  title:  “Approximation  of  multiparticle  Green’s  functions  via  sepa¬ 
rated  expansions”,  G.  Beylkin,  M.  Mohlenkamp,  L.  Monzon  and  F.  Perez 

(A  novel  approach  and  algorithms  for  constructing  accurate  approximations 
for  multiparticle  Green’s  functions  (for  energies  below  scattering  states)  using 
expansions  via  Gaussians) 

•  Tentative  title:  “Multiparticle  bound  states  of  Bose-Einstein  Condensate”,  G. 
Beylkin,  M.  Mohlenkamp,  L.  Monzon  and  F.  Perez 

(A  constructive  approach  to  computing  multiparticle  bound  states  for  confin¬ 
ing  quantum  harmonic  potential). 

•  Tentative  title:  “On  approximation  of  multi-variable  functions  by  exponential 
sums”,  G.  Beylkin  and  L.  Monzon 

(A  nonlinear  approximation  of  functions  of  several  variables  via  linear  com¬ 
bination  of  exponentials,  an  important  extension  of  [14]). 

•  Tentative  title:  “Approximating  a  wavefunction  as  an  unconstrained  sum  of 
Slater  determinants”,  G.  Beylkin,  M.  Mohlenkamp  and  F.  Perez 
(Methodology  for  using  [11]  for  quantum-mechanical  systems  and  basic  tests 
of  its  efficiency  and  accuracy.) 

•  Tentative  title:  “A  center-of-mass  principle  and  algorithmic  size-consistency 
for  the  multiparticle  Schrodinger  equation”,  G.  Beylkin  and  M.  Mohlenkamp 
(Structures  and  algorithms  for  size-consistency,  and  the  multipole-like  expan¬ 
sions  that  they  reveal.) 

•  Tentative  title:  “Capturing  the  interparticle  cusp  in  the  multiparticle  Schrodinger 
equation”,  G.  Beylkin  and  M.  Mohlenkamp 

(Structures  and  algorithms  for  simultaneously  capturing  the  cusps  between 
all  pairs  of  electrons.) 
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